Constant coefficients for the equation - Comparison with Central Difference scheme with IC and BCs.¶
In [145]:
from IPython.core.display import HTML
css_file = './custom.css'
HTML(open(css_file, "r").read())
Out[145]:
$$ \frac {\partial U}{\partial t} = D\frac {\partial^2 U}{\partial x^2} - v\frac {\partial U}{\partial x} + c*U + d $$$$ IC: U(x, 0) = \exp(-x**2) $$
$$ BC: Periodic: U(0, t) = U(1, t) $$ $$ \frac {\partial U(0, t)}{\partial x} = \frac {\partial U(1, t)}{\partial t} $$
In [130]:
from scipy.linalg import solve
# Define the equation parameters
from scipy.special import erfc
import pandas as pd
import numpy as np
import scipy.integrate
from numpy import exp
import matplotlib.pyplot as plt
In [ ]:
In [ ]:
In [131]:
def testAB_PS(a, b, c, d, nx, nt):
L = np.pi
T = 4.0
x = np.linspace(-L, L, nx+1)
t = np.linspace(0, T, nt+1)
u = np.exp(-x**2)#np.sin(np.pi*xin)
ui = u;
U_CDI = np.zeros((nx+1, nt+1))
error_2norm_CDI_PSTD = np.zeros(nt+1)
# Discretization of space and time
dx = 2*L/(nx+1)
# dt = 0.1*dx**2/D#T/(nt-1)
# dx = L/(nx+1)
dt = T/(nt+1)
#plt.figure()
#plt.plot(u, label='Initial')
def numerical_scheme_CDI(u_CDI, a, b, c, d, dx, dt):
alpha = a*dt/dx
beta = b*dt/dx**2
nx = len(u_CDI)-1;
A = np.zeros((nx+1, nx+1))
bi = d*dt+u_CDI;
for j in range(1, len(u_CDI)-1):
A[j, j-1] = -(beta-alpha/2)
A[j, j] = 1+2*beta-c*dt
A[j, j+1] = -(beta+alpha/2)
# Applying first order Periodic Boundary Condition
A[0, 0]=(1+2*beta);
A[0, 1]= (-alpha/2-beta)
A[0, nx]= (alpha/2-beta)
A[nx, 0]= (-alpha/2-beta) ;
A[nx, nx-1]= (alpha/2 - beta)
A[nx, nx]= 1+2*beta
# Applying second order periodic boundary conditions
# bt = d*dt+u_CDI;
# kk = np.zeros(2)
# bi = [*kk, *bt]
# A[0, 0]=0;
# A[0, 1]= 1
# A[0, N]= -1
# A[0, N+1]=0
# A[1, 0] = -1
# A[1,1] = 0
# A[1, 2] = 1
# A[1, N-1] = 1
# A[1, N] = 0
# A[1, N+1] = -1
# Boundary condition
# A[0, 0] = -3
# A[0, 1] = 4
# A[0, 2] = -1
# bi[0] = 0
# A[N, N] = 3
# A[N, N-1] = -4
# A[N, N-2] = 1
# bi[N] = 0
# A[0,0] = 1.0
# # A[N,N] = 1.0
# bi[0] = l
# bi[N] = r
u_CDI = solve(A, bi)
return u_CDI
# Apply the numerical scheme
# u = np.exp(-((x+0.5)**2)/(0.00125))
U_CDI = np.zeros((nx+1, nt+1))
error_2norm_QI_PSTD = np.zeros(nt)
u_CDI=u;
plt.plot(x,u)
for n in range(nt):
u_CDI = numerical_scheme_CDI(u_CDI, a, b, c, d, dx, dt)
U_CDI[:, n] = u_CDI
#plt.plot(u_CDI, label='Final')
#plt.show()
#plt.close()
plt.rcParams['font.size'] = 18
fig = plt.figure()
plt.imshow(U_CDI, cmap='viridis', extent=[0, T, 0, 2*L], aspect='auto')
plt.colorbar()
title = "CDI - nx = "+str(nx+1)
plt.title(title)
plt.xlabel('Time')
plt.ylabel('Position')
plt.show()
k = 2*np.pi*np.fft.rfftfreq(nx+1, 2*L/(nx))
k2=k**2;
# Defining variables
U_ps = np.zeros((nx+1, nt+1))
U_psk = np.zeros((nx+1, nt+1))
#plt.figure()
#plt.plot(xin, u, color = 'b', label = "initial")
# left = (0.025/(np.sqrt(0.000625 + 0.02*t)))*(np.exp((-(0.5 - t)**2)/(0.00125 + 0.04*t)))
# right = (0.025/(np.sqrt(0.000625 + 0.02*t)))*(np.exp((-(1.5 - t)**2)/(0.00125 + 0.04*t)))
# Solving over time
for i in range(nt+1): #
uk = np.fft.rfft(u)
uk[:] = (uk[:])/(1-(-b*k**2 + a*1j*k+c)*dt)
# print(np.shape(uk))
u = np.fft.irfft(uk) #+ d*dt
U_ps[:, i] = u
error = U_CDI[:, i] - u
error_2norm_CDI_PSTD[i] =np.linalg.norm(error)/np.sqrt(nx)
fig = plt.figure()
plt.imshow((U_ps), cmap='viridis', extent=[0, T, 0, 2*L], aspect='auto')
plt.colorbar()
plt.title("CD Scheme with PSTD")
plt.xlabel('Time')
plt.ylabel('Position')
title = "PSTDI - nx = "+str(nx+1)
plt.title(title)
plt.grid(True)
plt.show()
plt.close()
U_psk = np.zeros((nx+1, nt+1))
error_2norm_PSTDI_PSTDK = np.zeros(nt+1)
error_2norm_CDI_PSTDK = np.zeros(nt+1)
# Spatial grid
m=nx+1 # Number of grid points in space
L = np.pi # Width of spatial domain
x = np.arange(-m/2,m/2)*(2*L/m) # Grid points
dx = x[1]-x[0] # Grid spacing
# Temporal grid
tmax=T # Final time
N = nt+1 # number grid points in time
dt = tmax/N # interval between output times
xi = np.fft.fftfreq(m)*m*np.pi/L # Wavenumber "grid"
xi = 2*np.pi*np.fft.fftfreq(nx+1, 2*L/(nx))
# k = 2*np.pi*np.fft.rfftfreq(nx+1, 2*L/(nx))
# (this is the order in which numpy's FFT gives the frequencies)
# Initial data
u = np.exp(-x**2)# np.sin(2*x)**2 * (x<-L/4)
uhat0 = np.fft.fft(u)
plt.plot(x,u)
epsilon=b # Diffusion coefficient
# a = 1.0 # Advection coefficient
print(np.shape(xi))
print(np.shape(u))
# Now we solve the problem
for i in range(0,N):
t = i*dt
uhat = np.exp((1.j*xi*a - epsilon*xi**2)*t) * uhat0
u = np.real(np.fft.ifft(uhat)) #+ d*t
U_psk[:, i] = u
error = U_CDI[:, i] - u
error_2norm_CDI_PSTDK[i] =np.linalg.norm(error)/np.sqrt(nx)
error = U_ps[:, i] - u
error_2norm_PSTDI_PSTDK[i] =np.linalg.norm(error)/np.sqrt(nx)
fig = plt.figure()
plt.imshow((U_psk), cmap='viridis', extent=[0, T, 0, 2*L], aspect='auto')
plt.colorbar()
plt.title("CD Scheme with PSTDK")
plt.xlabel('Time')
plt.ylabel('Position')
title = "PSTDK - nx = "+str(nx+1)
plt.title(title)
plt.grid(True)
plt.show()
plt.close()
plt.figure();
# plt.plot(x, ui, '-', color='g', label = 'Initial');
plt.plot(x, U_CDI[:, int(nt/2)], '-*', color='b', label = 'CDI Method');
plt.plot(x, U_ps[:, int(nt/2)], '-+', color = 'r', label = 'Spectral Method');
plt.plot(x, U_psk[:, int(nt/2)], '-.', color = 'g', label = 'Ketchinson Spectral Method');
res = [U_CDI[:, nt-1], U_ps[:, nt-1], U_psk[:, nt-1]]
plt.legend(loc="upper right")
title = "Comparison - nx = "+str(nx+1)
plt.title(title)
plt.show()
plt.close()
plt.figure()
plt.plot(error_2norm_CDI_PSTD[0:nt], '.', label = 'PSTDI_CDI')
plt.plot(error_2norm_CDI_PSTDK[0:nt], '.', label = 'PSTDK_CDI')
plt.plot(error_2norm_PSTDI_PSTDK[0:nt], '.', label = 'PSTDI_PSTDK')
plt.legend()
plt.show()
plt.close()
return res
Pe = 100 without any source terms¶
Pe = 100¶
In [132]:
a = 1.0
b = 0.01
c = 0.0;
d = 0.0
testAB_PS(a, b, c, d, 63, 1000)
testAB_PS(a, b, c, d, 127, 1000)
res = testAB_PS(a, b, c, d, 255, 1000)
plt.figure();
plt.plot(res[0], '-*', color='b', label = 'CDI Method');
plt.plot(res[1], '-+', color = 'r', label = 'Spectral Method');
plt.plot(res[2], '-.', color = 'g', label = 'Ketchinson Spectral Method');
plt.legend(loc="upper right")
title = "Comparison - nx = "+str(255+1)
plt.title(title)
plt.show()
plt.close()
In [ ]:
Pe = 1¶
In [133]:
a = 1.0
b = 1.0
c = 0.0;
d = 0.0
testAB_PS(a, b, c, d, 63, 1000)
testAB_PS(a, b, c, d, 127, 1000)
testAB_PS(a, b, c, d, 255, 1000)
Out[133]:
Pe = 0.1¶
In [134]:
a = 0.1
b = 1.0
c = 0.0;
d = 0.0
testAB_PS(a, b, c, d, 63, 1000)
testAB_PS(a, b, c, d, 127, 1000)
testAB_PS(a, b, c, d, 255, 1000)
Out[134]:
In [135]:
a = 1.0
b = 1.0
c = 0.0
d = 0.0
print("sdeds")
testAB_PS(a, b, c, d, 63, 1000)
testAB_PS(a, b, c, d, 127, 1000)
testAB_PS(a, b, c, d, 255, 1000)
Out[135]:
In [136]:
a = 1.0
b = 0.01
c = 0.0;
d = 0.0
print("ssfw")
testAB_PS(a, b, c, d, 127, 1000)
testAB_PS(a, b, c, d, 127, 10000)
Out[136]:
With First Order Term¶
In [137]:
def testABC_PS(a, b, c, d, nx, nt):
L = np.pi
T = 4.0
x = np.linspace(-L, L, nx+1)
t = np.linspace(0, T, nt+1)
u = np.exp(-x**2)#np.sin(np.pi*xin)
ui = u;
U_CDI = np.zeros((nx+1, nt+1))
error_2norm_CDI_PSTD = np.zeros(nt+1)
# Discretization of space and time
dx = 2*L/(nx+1)
# dt = 0.1*dx**2/D#T/(nt-1)
# dx = L/(nx+1)
dt = T/(nt+1)
#plt.figure()
#plt.plot(u, label='Initial')
def numerical_scheme_CDI(u_CDI, a, b, c, d, dx, dt):
alpha = a*dt/dx
beta = b*dt/dx**2
nx = len(u_CDI)-1;
A = np.zeros((nx+1, nx+1))
bi = d*dt+u_CDI;
for j in range(1, len(u_CDI)-1):
A[j, j-1] = -(beta-alpha/2)
A[j, j] = 1+2*beta-c*dt
A[j, j+1] = -(beta+alpha/2)
# Applying first order Periodic Boundary Condition
A[0, 0]=(1+2*beta);
A[0, 1]= (-alpha/2-beta)
A[0, nx]= (alpha/2-beta)
A[nx, 0]= (-alpha/2-beta) ;
A[nx, nx-1]= (alpha/2 - beta)
A[nx, nx]= 1+2*beta
# Applying second order periodic boundary conditions
# bt = d*dt+u_CDI;
# kk = np.zeros(2)
# bi = [*kk, *bt]
# A[0, 0]=0;
# A[0, 1]= 1
# A[0, N]= -1
# A[0, N+1]=0
# A[1, 0] = -1
# A[1,1] = 0
# A[1, 2] = 1
# A[1, N-1] = 1
# A[1, N] = 0
# A[1, N+1] = -1
# Boundary condition
# A[0, 0] = -3
# A[0, 1] = 4
# A[0, 2] = -1
# bi[0] = 0
# A[N, N] = 3
# A[N, N-1] = -4
# A[N, N-2] = 1
# bi[N] = 0
# A[0,0] = 1.0
# # A[N,N] = 1.0
# bi[0] = l
# bi[N] = r
u_CDI = solve(A, bi)
return u_CDI
# Apply the numerical scheme
# u = np.exp(-((x+0.5)**2)/(0.00125))
U_CDI = np.zeros((nx+1, nt+1))
error_2norm_QI_PSTD = np.zeros(nt)
u_CDI=u;
plt.plot(x,u)
for n in range(nt):
u_CDI = numerical_scheme_CDI(u_CDI, a, b, c, d, dx, dt)
U_CDI[:, n] = u_CDI
#plt.plot(u_CDI, label='Final')
#plt.show()
#plt.close()
plt.rcParams['font.size'] = 18
fig = plt.figure()
plt.imshow(U_CDI, cmap='viridis', extent=[0, T, 0, 2*L], aspect='auto')
plt.colorbar()
title = "CDI - nx = "+str(nx+1)
plt.title(title)
plt.xlabel('Time')
plt.ylabel('Position')
plt.show()
k = 2*np.pi*np.fft.rfftfreq(nx+1, 2*L/(nx))
k2=k**2;
# Defining variables
U_ps = np.zeros((nx+1, nt+1))
U_psk = np.zeros((nx+1, nt+1))
#plt.figure()
#plt.plot(xin, u, color = 'b', label = "initial")
# left = (0.025/(np.sqrt(0.000625 + 0.02*t)))*(np.exp((-(0.5 - t)**2)/(0.00125 + 0.04*t)))
# right = (0.025/(np.sqrt(0.000625 + 0.02*t)))*(np.exp((-(1.5 - t)**2)/(0.00125 + 0.04*t)))
# Solving over time
for i in range(nt+1): #
uk = np.fft.rfft(u)
uk[:] = (uk[:])/(1-(-b*k**2 + a*1j*k+c)*dt)
# print(np.shape(uk))
u = np.fft.irfft(uk) #+ d*dt
U_ps[:, i] = u
error = U_CDI[:, i] - u
error_2norm_CDI_PSTD[i] =np.linalg.norm(error)/np.sqrt(nx)
fig = plt.figure()
plt.imshow((U_ps), cmap='viridis', extent=[0, T, 0, 2*L], aspect='auto')
plt.colorbar()
plt.title("CD Scheme with PSTD")
plt.xlabel('Time')
plt.ylabel('Position')
title = "PSTDI - nx = "+str(nx+1)
plt.title(title)
plt.grid(True)
plt.show()
plt.close()
U_psk = np.zeros((nx+1, nt+1))
error_2norm_PSTDI_PSTDK = np.zeros(nt+1)
error_2norm_CDI_PSTDK = np.zeros(nt+1)
# Spatial grid
m=nx+1 # Number of grid points in space
L = np.pi # Width of spatial domain
x = np.arange(-m/2,m/2)*(2*L/m) # Grid points
dx = x[1]-x[0] # Grid spacing
# Temporal grid
tmax=T # Final time
N = nt+1 # number grid points in time
dt = tmax/N # interval between output times
xi = np.fft.fftfreq(m)*m*np.pi/L # Wavenumber "grid"
xi = 2*np.pi*np.fft.fftfreq(nx+1, 2*L/(nx))
# k = 2*np.pi*np.fft.rfftfreq(nx+1, 2*L/(nx))
# (this is the order in which numpy's FFT gives the frequencies)
# Initial data
u = np.exp(-x**2)# np.sin(2*x)**2 * (x<-L/4)
uhat0 = np.fft.fft(u)
plt.plot(x,u)
epsilon=b # Diffusion coefficient
# a = 1.0 # Advection coefficient
print(np.shape(xi))
print(np.shape(u))
# Now we solve the problem
for i in range(0,N):
t = i*dt
uhat = np.exp((1.j*xi*a - epsilon*xi**2+c)*t) * uhat0
u = np.real(np.fft.ifft(uhat)) #+ d*t
U_psk[:, i] = u
error = U_CDI[:, i] - u
error_2norm_CDI_PSTDK[i] =np.linalg.norm(error)/np.sqrt(nx)
error = U_ps[:, i] - u
error_2norm_PSTDI_PSTDK[i] =np.linalg.norm(error)/np.sqrt(nx)
fig = plt.figure()
plt.imshow((U_psk), cmap='viridis', extent=[0, T, 0, 2*L], aspect='auto')
plt.colorbar()
plt.title("CD Scheme with PSTDK")
plt.xlabel('Time')
plt.ylabel('Position')
title = "PSTDK - nx = "+str(nx+1)
plt.title(title)
plt.grid(True)
plt.show()
plt.close()
plt.figure();
# plt.plot(x, ui, '-', color='g', label = 'Initial');
plt.plot(x, U_CDI[:, int(nt/2)], '-', color='b', label = 'CDI Method');
plt.plot(x, U_ps[:, int(nt/2)], '-', color = 'r', label = 'Spectral Method');
plt.plot(x, U_psk[:, int(nt/2)+1], '-', color = 'g', label = 'Ketchinson Spectral Method');
res = [U_CDI[:, nt-1], U_ps[:, nt-1], U_psk[:, nt]]
plt.legend(loc="upper right")
title = "Comparison - nx = "+str(nx+1)
plt.title(title)
plt.show()
plt.close()
plt.figure()
plt.plot(error_2norm_CDI_PSTD[0:nt], '.', label = 'PSTDI_CDI')
plt.plot(error_2norm_CDI_PSTDK[0:nt], '.', label = 'PSTDK_CDI')
plt.plot(error_2norm_PSTDI_PSTDK[0:nt], '.', label = 'PSTDI_PSTDK')
plt.legend()
plt.show()
plt.close()
return res
In [138]:
a = 1.0
b = 0.01
c = 0.5;
d = 0.0
testABC_PS(a, b, c, d, 63, 1000)
testABC_PS(a, b, c, d, 127, 1000)
res = testABC_PS(a, b, c, d, 255, 1000)
plt.figure();
plt.plot(res[0], '-', color='b', label = 'CDI Method');
plt.plot(res[1], '-', color = 'r', label = 'Spectral Method');
plt.plot(res[2], '-', color = 'g', label = 'Ketchinson Spectral Method');
plt.legend(loc="upper right")
title = "Comparison - nx = "+str(255+1)
plt.title(title)
plt.show()
plt.close()
In [139]:
a = 1.0
b = 1.0
c = 0.5;
d = 0.0
testABC_PS(a, b, c, d, 63, 1000)
testABC_PS(a, b, c, d, 127, 1000)
res = testABC_PS(a, b, c, d, 255, 1000)
plt.figure();
plt.plot(res[0], '-', color='b', label = 'CDI Method');
plt.plot(res[1], '-', color = 'r', label = 'Spectral Method');
plt.plot(res[2], '-', color = 'g', label = 'Ketchinson Spectral Method');
plt.legend(loc="upper right")
title = "Comparison - nx = "+str(255+1)
plt.title(title)
plt.show()
plt.close()
In [140]:
a = 0.1
b = 1.0
c = 0.5;
d = 0.0
testABC_PS(a, b, c, d, 63, 1000)
testABC_PS(a, b, c, d, 127, 1000)
res = testABC_PS(a, b, c, d, 255, 1000)
plt.figure();
plt.plot(res[0], '-', color='b', label = 'CDI Method');
plt.plot(res[1], '-', color = 'r', label = 'Spectral Method');
plt.plot(res[2], '-', color = 'g', label = 'Ketchinson Spectral Method');
plt.legend(loc="upper right")
title = "Comparison - nx = "+str(255+1)
plt.title(title)
plt.show()
plt.close()
With Zeroth and First Order Terms¶
In [141]:
def test_PS(a, b, c, d, nx, nt):
L = np.pi
T = 4.0
x = np.linspace(-L, L, nx+1)
t = np.linspace(0, T, nt+1)
u = np.exp(-x**2)#np.sin(np.pi*xin)
ui = u;
U_CDI = np.zeros((nx+1, nt+1))
error_2norm_CDI_PSTD = np.zeros(nt+1)
# Discretization of space and time
dx = 2*L/(nx+1)
# dt = 0.1*dx**2/D#T/(nt-1)
# dx = L/(nx+1)
dt = T/(nt+1)
#plt.figure()
#plt.plot(u, label='Initial')
def numerical_scheme_CDI(u_CDI, a, b, c, d, dx, dt):
alpha = a*dt/dx
beta = b*dt/dx**2
nx = len(u_CDI)-1;
A = np.zeros((nx+1, nx+1))
bi = d*dt+u_CDI;
for j in range(1, len(u_CDI)-1):
A[j, j-1] = -(beta-alpha/2)
A[j, j] = 1+2*beta-c*dt
A[j, j+1] = -(beta+alpha/2)
# Applying first order Periodic Boundary Condition
A[0, 0]=(1+2*beta);
A[0, 1]= (-alpha/2-beta)
A[0, nx]= (alpha/2-beta)
A[nx, 0]= (-alpha/2-beta) ;
A[nx, nx-1]= (alpha/2 - beta)
A[nx, nx]= 1+2*beta
# Applying second order periodic boundary conditions
# bt = d*dt+u_CDI;
# kk = np.zeros(2)
# bi = [*kk, *bt]
# A[0, 0]=0;
# A[0, 1]= 1
# A[0, N]= -1
# A[0, N+1]=0
# A[1, 0] = -1
# A[1,1] = 0
# A[1, 2] = 1
# A[1, N-1] = 1
# A[1, N] = 0
# A[1, N+1] = -1
# Boundary condition
# A[0, 0] = -3
# A[0, 1] = 4
# A[0, 2] = -1
# bi[0] = 0
# A[N, N] = 3
# A[N, N-1] = -4
# A[N, N-2] = 1
# bi[N] = 0
# A[0,0] = 1.0
# # A[N,N] = 1.0
# bi[0] = l
# bi[N] = r
u_CDI = solve(A, bi)
return u_CDI
# Apply the numerical scheme
# u = np.exp(-((x+0.5)**2)/(0.00125))
U_CDI = np.zeros((nx+1, nt+1))
error_2norm_QI_PSTD = np.zeros(nt)
u_CDI=u;
plt.plot(x,u)
for n in range(nt):
u_CDI = numerical_scheme_CDI(u_CDI, a, b, c, d, dx, dt)
U_CDI[:, n] = u_CDI
#plt.plot(u_CDI, label='Final')
#plt.show()
#plt.close()
plt.rcParams['font.size'] = 18
fig = plt.figure()
plt.imshow(U_CDI, cmap='viridis', extent=[0, T, 0, 2*L], aspect='auto')
plt.colorbar()
title = "CDI - nx = "+str(nx+1)
plt.title(title)
plt.xlabel('Time')
plt.ylabel('Position')
plt.show()
k = 2*np.pi*np.fft.rfftfreq(nx+1, 2*L/(nx))
k2=k**2;
# Defining variables
U_ps = np.zeros((nx+1, nt+1))
U_psk = np.zeros((nx+1, nt+1))
#plt.figure()
#plt.plot(xin, u, color = 'b', label = "initial")
# left = (0.025/(np.sqrt(0.000625 + 0.02*t)))*(np.exp((-(0.5 - t)**2)/(0.00125 + 0.04*t)))
# right = (0.025/(np.sqrt(0.000625 + 0.02*t)))*(np.exp((-(1.5 - t)**2)/(0.00125 + 0.04*t)))
# Solving over time
for i in range(nt+1): #
uk = np.fft.rfft(u)
uk[:] = (uk[:])/(1-(-b*k**2 + a*1j*k+c)*dt)
# print(np.shape(uk))
u = np.fft.irfft(uk) + d*dt
U_ps[:, i] = u
error = U_CDI[:, i] - u
error_2norm_CDI_PSTD[i] =np.linalg.norm(error)/np.sqrt(nx)
fig = plt.figure()
plt.imshow((U_ps), cmap='viridis', extent=[0, T, 0, 2*L], aspect='auto')
plt.colorbar()
plt.title("CD Scheme with PSTD")
plt.xlabel('Time')
plt.ylabel('Position')
title = "PSTDI - nx = "+str(nx+1)
plt.title(title)
plt.grid(True)
plt.show()
plt.close()
U_psk = np.zeros((nx+1, nt+1))
error_2norm_PSTDI_PSTDK = np.zeros(nt+1)
error_2norm_CDI_PSTDK = np.zeros(nt+1)
# Spatial grid
m=nx+1 # Number of grid points in space
L = np.pi # Width of spatial domain
x = np.arange(-m/2,m/2)*(2*L/m) # Grid points
dx = x[1]-x[0] # Grid spacing
# Temporal grid
tmax=T # Final time
N = nt+1 # number grid points in time
dt = tmax/N # interval between output times
xi = np.fft.fftfreq(m)*m*np.pi/L # Wavenumber "grid"
xi = 2*np.pi*np.fft.fftfreq(nx+1, 2*L/(nx))
# k = 2*np.pi*np.fft.rfftfreq(nx+1, 2*L/(nx))
# (this is the order in which numpy's FFT gives the frequencies)
# Initial data
u = np.exp(-x**2)# np.sin(2*x)**2 * (x<-L/4)
uhat0 = np.fft.fft(u)
plt.plot(x,u)
epsilon=b # Diffusion coefficient
# a = 1.0 # Advection coefficient
# print(np.shape(xi))
# print(np.shape(u))
# Now we solve the problem
d = np.ones(len(xi))*d
dfft = np.fft.fft(d)
# print(dfft)
# print((dfft/((1.j*xi*a - epsilon*xi**2+c)+1.j*xi)))
for i in range(0,N):
t = i*dt
al = uhat0 + dfft/((1.j*xi*a - epsilon*xi**2+c)+1.j*xi)
uhat = np.exp((1.j*xi*a - epsilon*xi**2+c)*t) * al-(dfft/((1.j*xi*a - epsilon*xi**2+c)+1.j*xi))
u = np.real(np.fft.ifft(uhat))
U_psk[:, i] = u
error = U_CDI[:, i] - u
error_2norm_CDI_PSTDK[i] =np.linalg.norm(error)/np.sqrt(nx)
error = U_ps[:, i] - u
error_2norm_PSTDI_PSTDK[i] =np.linalg.norm(error)/np.sqrt(nx)
plt.figure()
plt.imshow((U_psk), cmap='viridis', extent=[0, T, 0, 2*L], aspect='auto')
plt.colorbar()
plt.title("CD Scheme with PSTDK")
plt.xlabel('Time')
plt.ylabel('Position')
title = "PSTDK - nx = "+str(nx+1)
plt.title(title)
plt.grid(True)
plt.show()
plt.close()
plt.figure();
# plt.plot(x, ui, '-', color='g', label = 'Initial');
plt.plot(x, U_CDI[:, int(nt/2)], '-', color='b', label = 'CDI');
plt.plot(x, U_ps[:, int(nt/2)], '-', color = 'r', label = 'PSI');
plt.plot(x, U_psk[:, int(nt/2)+1], '-', color = 'g', label = 'KPS');
res = [U_CDI[:, nt-1], U_ps[:, nt-1], U_psk[:, nt]]
plt.legend(loc="upper right")
title = "Comparison - nx = "+str(nx+1)
plt.title(title)
plt.show()
plt.close()
plt.figure()
plt.plot(error_2norm_CDI_PSTD[0:nt], '.', label = 'PSTDI_CDI')
plt.plot(error_2norm_CDI_PSTDK[0:nt], '.', label = 'PSTDK_CDI')
plt.plot(error_2norm_PSTDI_PSTDK[0:nt], '.', label = 'PSTDI_PSTDK')
plt.legend()
plt.show()
plt.close()
return res
In [142]:
a = 1.0
b = 0.01
c = 0.1;
d = 0.5
test_PS(a, b, c, d, 63, 1000)
test_PS(a, b, c, d, 127, 1000)
res = test_PS(a, b, c, d, 255, 1000)
plt.figure();
plt.plot(res[0], '-', color='b', label = 'CDI Method');
plt.plot(res[1], '-', color = 'r', label = 'Spectral Method');
plt.plot(res[2], '-', color = 'g', label = 'Ketchinson Spectral Method');
plt.legend(loc="upper right")
title = "Comparison - nx = "+str(255+1)
plt.title(title)
plt.show()
plt.close()
In [143]:
a = 1.0
b = 1.0
c = 0.1;
d = 0.5
test_PS(a, b, c, d, 63, 1000)
test_PS(a, b, c, d, 127, 1000)
res = test_PS(a, b, c, d, 255, 1000)
plt.figure();
plt.plot(res[0], '-', color='b', label = 'CDI Method');
plt.plot(res[1], '-', color = 'r', label = 'Spectral Method');
plt.plot(res[2], '-', color = 'g', label = 'Ketchinson Spectral Method');
plt.legend(loc="upper right")
title = "Comparison - nx = "+str(255+1)
plt.title(title)
plt.show()
plt.close()
In [144]:
a = 0.1
b = 1.0
c = 0.1;
d = 0.5
test_PS(a, b, c, d, 63, 1000)
test_PS(a, b, c, d, 127, 1000)
res = test_PS(a, b, c, d, 255, 1000)
plt.figure();
plt.plot(res[0], '-', color='b', label = 'CDI Method');
plt.plot(res[1], '-', color = 'r', label = 'Spectral Method');
plt.plot(res[2], '-', color = 'g', label = 'Ketchinson Spectral Method');
plt.legend(loc="upper right")
title = "Comparison - nx = "+str(255+1)
plt.title(title)
plt.show()
plt.close()
Trying to compare with the semi-infinite solutions of Ogata and Banks¶
Solution for No Production Decay¶
Lapidus and Amundson 1952
Ogata and Banks 1961
Equation: R*dC/dt = D*d2C/dx2 - v*dC/dx
IC: C(x, 0) = Ci
BCs:Dirichlet for 0<t<to: C(0, t) = Co
Neumann: dC/dx(inf, t) = 0
$$ BC - Left: U(0, t) = C_o \ \ \ \ 0<t<t_o $$ $$ BC - Right: \frac {dU(\infty, t)}{dx} = 0 $$ $$ c = 0 $$ $$ d = 0 $$
In [26]:
def Analytical_Ogata(a, b, c, d, L, T, nx, nt):
# Analytical - du/dt = -v*du/dx+D*d2u/dx2
Ci = 0;
Co = 2;
# dCinf/dx = 0;
R = 1;
D = b
v = -a
print(v, D, c, d)
x = np.linspace(0, L, nx+1);
t = np.linspace(0, T, nt);
C = np.zeros((nt, nx))
for i in range(nt):
for j in range(nx):
C[i, j] = Ci + ((Co-Ci)/2)*(erfc((x[j]-v*t[i])/(2*np.sqrt(D*t[i]/R))) + (np.exp(v*x[j]/(D/R)))*(erfc((x[j]+v*t[i])/(2*np.sqrt(D*t[i]/R)))))
U_exact = np.transpose(C);
fig = plt.figure()
plt.imshow(np.transpose(C), cmap='viridis', extent=[0, T, L, 0], aspect='auto')
plt.colorbar()
plt.title("Exact")
plt.xlabel('Time')
plt.ylabel('Position')
#fig.set_size_inches(30.,18.)
# plt.savefig('U_CDI.png', dpi = 900)
plt.show()
return U_exact[:, -1]
In [27]:
## Ogata and Banks - No production or decay terms
axs = [-1.0, -0.1, -0.01]
bs = [0.01, 0.1]
it = 0
L = 1.0;
T = 5.0;
nx = 128;
nt = 5001;
plt.figure()
print("ddfed")
U_compab = np.zeros((12, nx))
for a in axs:
for b in bs:
pe = abs(a/b);
# U_exact =
U_compab[it, :] = Analytical_Ogata(a, b, 0, 0, L, T, nx, nt)
U_compab[it+1, :] = testCD_PS(-a, b, 0, 0, nx-1, nt-1)
# U_comp = comparePeclet(a, b)
plt.plot(U_compab[i, :],'+-'); plt.plot(U_compab[i+1, :],'-', label = str(np.round(pe, 3)));
plt.ylim([0, 2.1])
# plt.title("Pe = "+str(np.round(pe, 3)))
plt.xlabel("Dimensionless Time")
plt.ylabel("Dimensionless Concentration")
plt.title("Ogata and Banks - Different Pe="+str(np.round(pe, 3)))
it = it+2
# Printing the temporal grid tests
zipped1 = list(zip(U_compab[0, :], U_compab[1, :], U_compab[2, :], U_compab[3, :], U_compab[4, :], U_compab[5, :], U_compab[6, :], U_compab[7, :], U_compab[8, :], U_compab[9, :], U_compab[10, :], U_compab[11, :]))
df2 = pd.DataFrame(zipped1, columns=['OB', 'Pe = 100.0', 'OB', 'Pe = 10.0', 'OB', 'Pe = 10.0','OB', 'Pe = 1.0','OB', 'Pe = 1.0','OB', 'Pe = 0.1'])
df2.to_csv('Analysis_OgataBanks_complete.csv', index=False)
Didn't work - Using MMS¶
MMS for the same manufactured solution as it satisfies the present condition¶
For an equation with given IC and periodic BCs : $$ \frac {\partial U}{\partial t} = D\frac {\partial^2 U}{\partial x^2} + v\frac {\partial U}{\partial x}$$ $$ IC: U(x, 0) = C_i $$ $$ BC - Left: U(0, t) = -C_o $$ $$ BC - Right: U(0, t) = -C_o $$ $$ BC - Left: \frac {dU(0, t)}{dx} = 0 $$ $$ BC - Right: \frac {dU(1, t)}{dx} = 0 $$
The chosen equation is: $$ U(x,t) = C_o \cos(\pi(2x -1)) e^{-t}$$
The initial and boundary condition terms: I.C.: Cisin(pix) Boundary Conditions: Periodic BC $$ U(x, 0) = C_o\cos( \pi(2x -1))$$ $$ U(0, t) = -C_o e^{-t} $$ $$ U(1, t) = -C_o e^{-t} $$
$$ \frac {dU}{dx} (0, t) = -2\pi C_o\sin( \pi(2*0 -1)) e^{-t} = 0$$$$ \frac {dU}{dx} (1, t) = -2\pi C_o\sin(\pi(2*1 -1)) e^{-t} = 0$$
The derivative terms are:
$$ \frac {dU}{dt} = -C_o\cos(\pi (2*x-1)) e^{-t}$$$$ \frac {dU}{dx} = -C_o 2 \pi \sin(\pi (2*x-1)) e^{-t}$$$$ \frac {d^2U}{dx^2} = -C_o 4 \pi^2 \cos(\pi (2*x-1)) e^{-t}$$The new source term will be:
$$ d = \frac {dU}{dt} - v*\frac {dU}{dx} - D*\frac {d^2U}{dx^2}$$$$ d = -C_o\cos(\pi (2*x-1)) e^{-t} - v*2*C_o\pi \sin(\pi (2*x-1)) e^{-t} - D*C_o*4*\pi^2 \cos(\pi (2*x-1)) e^{-t}$$$$ d = -C_o e^{-t}(\cos(\pi (2*x-1)) - v*2*\pi \sin(\pi (2*x-1)) - D*4*\pi^2 \cos(\pi (2*x-1))) $$In [28]:
def numerical_scheme_CDI(u_CDI, a, b, c, d, dx, dt):
alpha = b*dt/dx
beta = a*dt/dx**2
nx = len(u_CDI)-1;
A = np.zeros((nx+1, nx+1))
bi = d*dt+u_CDI;
for j in range(1, len(u_CDI)-1):
A[j, j-1] = -(beta-alpha/2)
A[j, j] = 1+2*beta-c*dt
A[j, j+1] = -(beta+alpha/2)
# Applying first order Periodic Boundary Condition
A[0, 0]=(1+2*beta)-c*dt;
A[0, 1]= (-alpha/2-beta)
A[0, nx]= (alpha/2-beta)
A[nx, 0]= (-alpha/2-beta) ;
A[nx, nx-1]= (alpha/2 - beta)
A[nx, nx]= 1+2*beta-c*dt
u_CDI = solve(A, bi)
return u_CDI
In [29]:
def numerical_scheme_PSTD(u, a, b, c, d, k, dx, dt):
uk = np.fft.rfft(u)
uk[:] = (uk[:])/(1-(-a*k**2 + b*1j*k+c)*dt)
u = np.fft.irfft(uk) + d*dt
return u
In [30]:
# Creating function for running the code and getting the error with different number of nodes with left dirichlet and right neumann condition.
def numericalcodes(nx, nt):
# Define the grid
a = 1.0;
b = 0.01;
c = 0.0;
Co = 1.0;
L = 1.0
T = 1.0
dm = np.zeros((nx+1, nt+1));
U_exact = np.zeros((nx+1, nt+1));
# Define x & t
x = np.linspace(0, L, nx+1)
t = np.linspace(0, T, nt+1)
for i in range(nt+1):
for j in range (nx+1):
dm[j, i] = -(Co*np.exp(-t[i]))*(np.cos(np.pi*(2*x[j]-1)) - 2*b*np.pi*np.sin(np.pi*(2*x[j]-1)) - a*4*(np.pi**2)*np.cos(np.pi*(2*x[j]-1)))
for i in range(nt+1):
for j in range (nx+1):
U_exact[j, i] = Co*np.cos(np.pi*(2*x[j]-1))*np.exp(-t[i])
fig = plt.figure()
plt.imshow((U_exact), cmap='viridis', extent=[0, L, 0, T], aspect='auto')
plt.colorbar()
plt.title("Exact")
plt.xlabel('Time')
plt.ylabel('Position')
plt.show()
# Defining variables and errors
U_CDI = np.zeros((nx+1, nt+1))
U_ps = np.zeros((nx+1, nt+1))
error_2norm_CDI_PSTD = np.zeros(nt+1)
error_2norm_Exact_CDI = np.zeros(nt+1)
error_2norm_Exact_PSTD = np.zeros(nt+1)
# Discretization of space and time
dx = L/(nx)
dt = T/(nt)
print("\n here \n")
print(dx)
print(x[1])
print(dt)
print(t[1])
print("\n here \n")
u = Co*np.cos(np.pi*(2*x-1)) # np.zeros(nx);
u_ps = u;
u_CDI=u;
u1 = u
k = 2*np.pi*np.fft.rfftfreq(nx+1, L/(nx))
# Solving over time
for i in range(nt+1):
u_CDI = numerical_scheme_CDI(u_CDI, a, b, c, dm[:, i], dx, dt)
u_ps = numerical_scheme_PSTD(u_ps, a, b, c, dm[:, i], k, dx, dt)
if(i<2):
u1= u_ps
U_CDI[:, i] = u_CDI
U_ps[:, i] = u_ps
error = U_CDI[:, i] - U_ps[:, i]
error_2norm_CDI_PSTD[i] =np.linalg.norm(error)/np.sqrt(nx)
error = U_ps[:, i] - U_exact[:, i]
error_2norm_Exact_PSTD[i] =np.linalg.norm(error)/np.sqrt(nx)
errorPSTD = error_2norm_Exact_PSTD[i]
error = U_CDI[:, i] - U_exact[:, i]
error_2norm_Exact_CDI[i] =np.linalg.norm(error)/np.sqrt(nx)
errorCDI = error_2norm_Exact_CDI[i]
# Errors at time 1s for the two methods(CDI, PSTD) and 7 number of nodes (dx in 20, 40, 60, 80, 100, 120, 150)
errornorms = np.ones(7)
errornorms[0] = errorCDI
errornorms[1] = errorPSTD
plt.rcParams['font.size'] = 18
fig = plt.figure()
plt.imshow(U_CDI, cmap='viridis', extent=[0, T, 0, L], aspect='auto')
plt.colorbar()
title = "CDI - nx = "+str(nx+1)
plt.title(title)
plt.xlabel('Time')
plt.ylabel('Position')
plt.show()
plt.figure()
plt.imshow((U_ps), cmap='viridis', extent=[0, T, 0, L], aspect='auto')
plt.colorbar()
plt.title("CD Scheme with PSTD")
plt.xlabel('Time')
plt.ylabel('Position')
title = "PSTDI - nx = "+str(nx+1)
plt.title(title)
plt.grid(True)
plt.show()
plt.close()
plt.figure();
plt.plot(x, u, '-', color='g', label = 'Initial');
plt.plot(x, u1, '-', color='c', label = 'InitialPS');
plt.plot(x, U_CDI[:, nt-1], '-*', color='b', label = 'CDI Method');
plt.plot(x, U_ps[:, nt-1], '-+', color = 'r', label = 'Spectral Method');
plt.plot(x, U_exact[:, nt-1], '.-', color = 'k', label = 'Exact');
plt.legend(loc="upper right")
title = "Comparison - nx = "+str(nx+1)
plt.title(title)
plt.show()
plt.close()
plt.figure()
plt.plot(error_2norm_CDI_PSTD[0:nt+1], '.', label = 'PSTDI_CDI')
plt.plot(error_2norm_Exact_PSTD[0:nt+1], '.', label = 'Exact_PSTD')
plt.plot(error_2norm_Exact_CDI[0:nt+1], '.', label = 'Exact_CDI')
plt.legend()
plt.show()
plt.close()
return errornorms
In [31]:
# Equation with decay terms
# Spatial MMS tests
errorsx = np.ones((8, 2))
errornorms = np.ones(8)
nx = 32;
nt = 1000;
nxs = [15, 31, 63, 127, 255, 511, 1023, 2047]
nts = [250, 500, 1000, 2000, 4000]
for i in range(6):
errornorms = numericalcodes(nxs[i], nt)
errorsx[i, 0] = errornorms[0]
errorsx[i, 1] = errornorms[1]
fig = plt.figure()
plt.figure()
fig.set_size_inches(10.,2.)
plt.plot(nxs, errorsx[:, 0], '-', label = 'CDI')
plt.plot(nxs, errorsx[:, 1], '.-', label = 'PSTD')
plt.xlabel('Grid Numbers')
plt.ylabel('Error Norms')
plt.legend()
fig = plt.figure()
fig.set_size_inches(10.,7.)
px = np.zeros((6, 2))
for g in range(6):
for k in range(2):
px[g, k] = np.log((np.abs(errorsx[g+2, k]) - np.abs(errorsx[g+1, k]))/(np.abs(errorsx[g+1, k]) - np.abs(errorsx[g, k])))/(np.log(nxs[g]/nxs[g+1]))
plt.plot(nxs[0:6], px[:, 0], '-', label = 'CDI')
plt.plot(nxs[0:6], px[:, 1], '-.', label = 'PSTD')
plt.xlabel('Spatial Grid Numbers')
plt.ylabel('Observed Order')
plt.legend()
# Printing the spatial grid tests
zipped1 = list(zip(nxs, errorsx[:, 0], errorsx[:, 1], px[:, 0], px[:, 1]))
df2 = pd.DataFrame(zipped1, columns=['Spatial Discretization Numbers', 'CDI_Error', 'PSTD_Error', 'CDI_Order', 'PSTD_Order'])
df2.to_csv('Analysis_Spatial25_Complete.csv', index=False)
In [32]:
a = 1.0;
b = 0.01;
c = 0.0;
Co = 1.0;
L = 1.0
T = 1.0
dm = np.zeros((nx+1, nt+1));
U_exact = np.zeros((nx+1, nt+1));
# Define x & t
x = np.linspace(0, L, nx+1)
t = np.linspace(0, T, nt+1)
for i in range(nt+1):
for j in range (nx+1):
U_exact[j, i] = Co*np.cos(np.pi*(2*x[j]-1))*np.exp(-t[i])
plt.plot(x, U_exact[:, 1], '.-', color = 'k', label = 'Exact0');
plt.plot(x, U_exact[:, 500], '.-', color = 'g', label = 'Exact500');
plt.plot(x, U_exact[:, 900], '.-', color = 'b', label = 'Exact900');
In [33]:
def tests2(a, b, c, d, nx, nt):
L = 1.0
T = 1.0
x = np.linspace(0, L, nx+1)
t = np.linspace(0, T, nt+1)
xin = np.linspace(0, L, nx+1);
u = np.exp(-x**2)#np.sin(np.pi*xin)
ui = u;
U_CDI = np.zeros((nx+1, nt+1))
error_2norm_CDI_PSTD = np.zeros(nt+1)
# Discretization of space and time
dx = L/(nx-1)
dt = T/(nt+1)
#plt.figure()
#plt.plot(u, label='Initial')
def numerical_scheme_CDI(u_CDI, a, b, c, d, dx, dt):
alpha = a*dt/dx
beta = b*dt/dx**2
nx = len(u_CDI)-1;
A = np.zeros((nx+1, nx+1))
bi = d*dt+u_CDI;
for j in range(1, len(u_CDI)-1):
A[j, j-1] = -(beta[j]-alpha[j]/2)
A[j, j] = 1+2*beta[j]-c[j]*dt
A[j, j+1] = -(beta[j]+alpha[j]/2)
# Applying first order Periodic Boundary Condition
A[0, 0]=(1+2*beta[0])-c[0]*dt;
A[0, 1]= (-alpha[0]/2-beta[0])
A[0, nx]= (alpha[0]/2-beta[0])
A[nx, 0]= (-alpha[nx]/2-beta[nx]) ;
A[nx, nx-1]= (alpha[nx]/2 - beta[nx])
A[nx, nx]= 1+2*beta[nx]-c[nx]*dt
u_CDI = solve(A, bi)
return u_CDI
# Apply the numerical scheme
U_CDI = np.zeros((nx+1, nt+1))
error_2norm_QI_PSTD = np.zeros(nt)
u_CDI=u;
for n in range(nt):
u_CDI = numerical_scheme_CDI(u_CDI, a, b, c, d, dx, dt)
U_CDI[:, n] = u_CDI
plt.rcParams['font.size'] = 18
fig = plt.figure()
plt.imshow(U_CDI, cmap='viridis', extent=[0, T, 0, L], aspect='auto')
plt.colorbar()
title = "CDI - nx = "+str(nx+1)
plt.title(title)
plt.xlabel('Time')
plt.ylabel('Position')
plt.show()
k = 2*np.pi*np.fft.rfftfreq(nx+1, L/(nx))
k2=k**2;
# Defining variables
U_ps = np.zeros((nx+1, nt+1))
# Solving over time
for i in range(nt+1): #
uk = np.fft.rfft(u)
uk[:] = (uk[:])/(1-(-b*k**2 + a*1j*k+c)*dt)
u = np.fft.irfft(uk) + d*dt
U_ps[:, i] = u
error = U_CDI[:, i] - u
error_2norm_CDI_PSTD[i] =np.linalg.norm(error)/np.sqrt(nx)
fig = plt.figure()
plt.imshow((U_ps), cmap='viridis', extent=[0, T, 0, L], aspect='auto')
plt.colorbar()
plt.title("CD Scheme with PSTD")
plt.xlabel('Time')
plt.ylabel('Position')
title = "PSTDI - nx = "+str(nx+1)
plt.title(title)
plt.grid(True)
plt.show()
plt.close()
plt.figure();
plt.plot(x, ui, '-', color='g', label = 'Initial');
plt.plot(x, U_CDI[:, nt-1], '-*', color='b', label = 'CDI Method');
plt.plot(x, U_ps[:, nt-1], '-+', color = 'r', label = 'Spectral Method');
res = [U_CDI[:, nt-1], U_ps[:, nt-1]]
# print(np.shape(res))
plt.legend(loc="upper right")
title = "Comparison - nx = "+str(nx+1)
plt.title(title)
plt.show()
plt.close()
plt.figure()
plt.plot(error_2norm_CDI_PSTD[0:nt], '.', label = 'PSTDI_CDI')
plt.legend()
plt.show()
plt.close()
return res
In [34]:
nx = 63
x = np.linspace(0,1, nx+1)
a = 1.0*np.sin(np.pi*x)
b = 0.01*np.cos(np.pi*x)
c = a+10*b;
d = a-10*b;
tests2(a, b, c, d, nx, 1000)
nx = 127
x = np.linspace(0,1, nx+1)
a = 1.0*np.sin(np.pi*x)
b = 0.01*np.cos(np.pi*x)
c = a+10*b;
d = a-10*b;
tests2(a, b, c, d, nx, 1000)
nx = 255
x = np.linspace(0,1, nx+1)
a = 1.0*np.sin(np.pi*x)
b = 0.01*np.cos(np.pi*x)
c = a+10*b;
d = a-10*b;
res = tests2(a, b, c, d, nx, 1000)
print('ldd')
plt.figure();
plt.plot(res[0], '-*', color='b', label = 'CDI Method');
plt.plot(res[1], '-+', color = 'r', label = 'Spectral Method');
plt.legend(loc="upper right")
title = "Comparison - nx = "+str(nx+1)
plt.title(title)
plt.show()
plt.close()
In [ ]:
In [ ]:
In [ ]:
In [ ]: